Method of calculating optical path distribution inside scattering absorber

ABSTRACT

A method comprises a first step of dividing a measured medium into N voxels; a second step of assuming an imaginary medium identical to the measured medium except for absence of absorption; a third step of setting a light injection position and a light detection position for the imaginary medium; a fourth step of calculating an impulse response s(t) at the light detection position; a fifth step of calculating a probability density Ui(t′,t) that the impulse response at a time t at the light detection position has existed in voxel i (i=1, . . . , N; 1≦N) at a point of time t′ (0≦t′≦t); a sixth step of calculating a cumulative probability Ui(t) that the impulse response detected at the time t has existed in voxel i, from Ui(t′,t); and an eighth step of calculating a time-resolved path length li in voxel i from Ui(t) and s(t).

TECHNICAL FIELD

[0001] The present invention relates to methods of calculating adistribution of photon paths inside a scattering medium. Moreparticularly, the invention relates to methods of calculating atime-resolved photon path distribution inside a scattering medium, whichcan be applied to apparatus configured to move a light injectionposition and a light detection position along a surface of a measuredobject and gain information about the interior of the measured object.

[0002] From this time-resolved photon path distribution, we cancalculate a distribution of photon paths in an appropriate measuringmethod and calculate contributions of respective portions of thescattering medium to an amount of optical attenuation or to a mean pathlength (also referred to as a contribution function or weight function).

BACKGROUND ART

[0003] In the field of optical CT or the like in a broad sense,including interior measurement, there are conventionally known measuringmethods of gaining information about the interior of the measuredobject, using the weight function or contribution function, and a numberof reports have already been made about the weight functions and othersapplied thereto. These are described, for example, in the followingdocuments.

[0004] (1) S. Arridge: SPIE Institutes for Advanced OpticalTechnologies, Vol. IS11, Medical Optical Tomography: Functional Imagingand Monitoring, 35-64(1993); (2) R. L. Barbour and H. L. Graber: ibid.87-120(1993); (3) H. L. Graber, J. Chang, R. Aronson and R. L. Barbour:ibid. 121-143(1993); (4) J. C. Schotland, J. C. Haselgrove and J. S.Leigh: Applied Optics, 32, 448-5453(1993); (5) Chang, R. Aronson, H. L.Graber and R. L. Barbour: Proc. SPIE, 2389, 448-464(1995); (6) B. W.Pogue, M. S. Patterson, H. Jiang and K. D. Paulsen: Phys. Med. Biol. 40,1709-1729(1995); (7) S. R. Arridge: Applied Optics, 34, 7395-7409(1995);(8) H. L. Graber, J. Chang, and R. L. Barbour: Proc. SPIE, 2570,219-234(1995); (9) A. Maki and H. Koizumi: OSA TOPS, Vol.2,299-304(1996); (10) H. Jiang, K. D. Paulsen and Ulf L. Osterberg: J.Opt. Soc. Am. A13, 253-266(1996); (11) S. R. Arridge and J. C. Hebden:Phys. Med. Biol. 42, 841-853(1997); (12) S. B. Colak, D. G. Papaioannou,G. W. 't Hooft, M. B. van der Mark, H. Schomberg, J. C. J. Paasschens,J. B. M. Melissen and N. A. A. J. van Astten: Applied Optics, 36,180-213(1997).

[0005] However, the interior information measuring methods of opticalCTs and others disclosed in the above documents (1) to (12) and theweight functions applied thereto included the problems discussed belowand thus posed significant problems in measurement accuracy andapplicability. In fact, no report has been made yet about an example ofpractical use of optical CT with satisfactory performance in theforegoing field.

[0006] Specifically, the first problem in the conventional optical CTand others is that the photon migration analysis or photon migrationmodel in the medium is based on the photon diffusion equation asapplication of diffusion approximation to the transport equation. Inother words, since the diffusion approximation holds only in mediasufficiently larger than the mean free pathlength of photons in themedium, there is a problem that it is impossible to handle relativelysmall media, tissue of complicated interior shape, and media ofcomplicated shape. In addition, the diffusion approximation is premisedon isotropic scattering, and thus application thereof to the measurementof the measured object like living tissue having anisotropic scatteringcharacteristics will pose a problem that there occur considerable errorsdue to the approximation of anisotropic scattering to isotropicscattering.

[0007] Further, differential equations such as the diffusion equationand the like involve a bothersome problem that even with any numericalcomputation approach such as the analytical or finite element method,boundary conditions (the shape of the medium, reflection characteristicsat interfaces, etc.) must be preliminarily set and then a solution canbe determined. Namely, in the case of the measured object like livingtissue, the boundary conditions normally vary depending upon a place tobe measured, the wavelength of light used in measurement, and so on, andfor improvement in accuracy on the basis of correction for influence ofthese factors, it is necessary to repeat complicated calculations atevery change of the boundary conditions, which results in a big problemof extremely long calculation time.

[0008] The second problem is that the weight function in a narrow sense(also called the contribution function), i.e., a mean path length(weighted mean) for an ensemble of photons constituting an impulseresponse of the medium, or a phase delay equivalent thereto (measured inthe frequency domain) is applied to the calculation of the informationon the interior of the measured object. In this case, the weightfunction in a narrow sense (the mean path length or the phase delayequivalent thereto) varies depending upon absorption coefficients andabsorption distribution, and thus handling thereof is extremelycomplicated. Since practical use of a calculation method taking accountof such dependence inevitably requires calculations in a large iterationnumber, there will arise a problem that the calculation time becomesextremely long over a practical range. Therefore, the dependence on theabsorption coefficients and absorption distribution is normally ignored,but this approximation can be the cause of a serious issue of increasein errors.

[0009] In order to reduce the errors due to the dependence of the weightfunction on the absorption coefficients and absorption distribution asdescribed above, there is a method of calculating and applying theweight function with some appropriate absorption. However, this poses aproblem that the time necessary for the calculation of the weightfunction with absorption is extremely longer than the time necessary forthe calculation of the weight function without absorption.

[0010] For the reason described above, it is concluded that theconventional measuring methods using the weight function in a narrowsense (the mean path length or the phase delay equivalent thereto) arenot practical.

[0011] Besides the application of the above weight function in a narrowsense, there is a further measuring method of applying the perturbationtheory to the approximate equation of the transport equation or to thephoton diffusion equation to gain the information about the interior ofthe measured object with the use of the relationship between signallight and optical characteristics of the scattering medium. However,this method requires extremely complex handling of the nonlinear effects(second and higher order terms). In this respect, the calculationincluding the second and higher order terms can be theoreticallyperformed by a computer, but the calculation time will be huge even withthe use of the currently fastest computer, so as to make practical usethereof impossible. It is thus common practice to ignore the second andhigher order terms. For this reason, this method had a problem thatthere arose large errors due to interaction between absorbing regions inreconstruction of an optical CT image of a medium containing a pluralityof relatively strong absorbing regions.

[0012] As described above, the measuring methods of gaining theinformation about the interior of the measured object, such as theconventional optical CT and the like, failed to obtain a reconstructedimage with satisfactory accuracy and had the significant issues as tothe spatial resolution, image distortion, quantitation, measurementsensitivity, required measurement time, and so on, which made practicaluse thereof difficult.

[0013] The Inventor contemplated that the following was important inorder to break through the above circumstances, and has pursued a seriesof studies. Namely, an important subject for achievement of a measuringmethod capable of obtaining the information about the interior of themeasured object, particularly, for achievement of optical CT, is toclarify the details of behavior of photons migrating in living tissue asa strong scattering medium, more precisely describe the relation betweendetected signal light and the optical characteristics of the scatteringmedium (scattering absorber) containing an absorptive constituent, anddevelop a new algorithm of reconstructing an optical CT image by makinguse of the signal light and the relation between the signal light andthe optical characteristics of the scattering medium. Then the Inventorcontemplated applying the Microscopic Beer-Lambert Law (hereinafterreferred to as “MBL”) to the analysis of the behavior of photonsmigrating in the scattering medium and making use of information aboutthe photon path distribution, i.e., about where photons have passed, forthe reconstruction of optical CT image.

[0014] Then the Inventor has proposed a photon migration model (a finitegrid model) corresponding to the scattering media, based on the MBL,derived an analytic equation representing the relation between theoptical characteristics of the scattering media and the signal light,and developed methods of analyzing the behavior of photons in thescattering media. The Inventor and others reported the results of these,for example, in the following documents.

[0015] (13) Y. Tsuchiya and T. Urakami: “Photon migration model forturbid biological medium having various shapes”, Jpn. J. Appl. Phys. 34.Part2, pp. L79-81(1995); (14) Y. Tsuchiya and T. Urakami, “Frequencydomain analysis of photon migration based on the microscopicBeer-Lambert law”, Jpn. J. Appl. Phys. 35, Part1, pp.4848-4851(1996);(15) Y. Tsuchiya and T. Urakami: “Non-invasive spectroscopy of variouslyshaped turbid media like human tissue based on the microscopicBeer-Lambert law”, OSA TOPS, Biomedical Optical Spectroscopy andDiagnostics 1996, 3, pp.98-100(1996); (16) Y. Tsuchiya and T. Urakami:“Quantitation of absorbing substances in turbid media such as humantissues based on the microscopic Beer-Lambert law”. Optics Commun. 144,pp.269-280(1997); (17) Y. Tsuchiya and T. Urakami: “Optical quantitationof absorbers in variously shaped turbid media based on the microscopicBeer-Lambert law: A new approach to optical computerized tomography”,Advances in Optical Biopsy and Optical Mammography (Annals of the NewYork Academy of Sciences), 838, pp.75-94(1998); (18) Y. Ueda, K. Ohta,M. Oda, M. Miwa, Y. Yamashita, and Y. Tsuchiya: “Average value method: Anew approach to practical optical computed tomography for a turbidmedium such as human tissue”, Jpn. J. Appl. Phys. 37, Part1. 5A,pp.2717-2723(1998); (19) Yutaka Tsuchiya; “Reconstruction of optical CTimage based on Microscopic Beer-Lambert Law and average value method,” Oplus E, Vol.21, No.7, 814-821; (20) H. Zhang, M. Miwa, Y. Yamashita, andY. Tsuchiya: “Quantitation of absorbers in turbid media usingtime-integrated spectroscopy based on microscopic Beer-Lambert law”,Jpn. J. Appl. Phys. 37, Part1, pp.2724-2727(1998); (21) H. Zhang, Y.Tsuchiya, T. Urakami, M. Miwa, and Y. Yamashita: “Time integratedspectroscopy of turbid media based on the microscopic Beer-Lambert law:Consideration of the wavelength dependence of scattering properties”,Optics Commun. 153, pp.314-322(1998); (22) Y. Tsuchiya, H. Zhang, T.Urakami, M. Miwa, and Y. Yamashita: “Time integrated spectroscopy ofturbid media such as human tissues based on the microscopic Beer-Lambertlaw”, Proc. JICAST'98/CPST'98, Joint international conference onAdvanced science and technology, Hamamatsu, August 29-30, pp.237-240(1998); (23) H. Zhang, T. Urakami, Y. Tsuchiya, Z. Liu, and T.Hiruma: “Time integrated spectroscopy of turbid media based on themicroscopic Beer-Lambert law: Application to small-size phantoms havingdifferent boundary conditions”, J. Biomedical Optics. 4,pp.183-190(1999);(24) Y. Tsuchiya, Y. Ueda, H. Zhang, Y. Yamashita, M.Oda, and T. Urakami: “Analytic expressions for determining theconcentrations of absorber in turbid media by time-gating measurements”,OSA TOPS, Advances in Optical Imaging and Photon Migration, 21,pp.67-72(1998); (25) H. Zhang, M. Miwa, T. Urakami, Y. Yamashita, and Y.Tsuchiya: “Simple subtraction method for determining the mean pathlength traveled by photons in turbid media”, Jpn. J. Appl. Phys. 37,Part1, pp.700-704 (1998).

[0016] This MBL is the law that “in a microscopic view in a mediumhaving arbitrary scattering and absorption distributions, a photonmigrating in a portion having the absorption coefficient of μ_(a), isexponentially attenuated because of absorption along a migrating zigzagphoton path of length l, a survival probability of the photon is givenby the value exp(−μ_(a)l) independent of the scattering characteristicsof the medium and the boundary conditions, and the amount of attenuationdue to the absorption is μ_(a)l,” and is expressed, for example, by thefollowing equation.

ƒ(l)=ƒ⁰(l)exp(−μ_(a)l)   (1)

[0017] Here f(l) is a detection probability density function of photonsactually measured, and f⁰(l) a detection probability density function ofphotons in a state of no absorption.

[0018] This Eq (1) indicates that in the medium having arbitraryscattering and absorption distributions, absorption and scatteringevents are independent of each other and that the superpositionprinciple holds as to the absorption, i.e., that, in the case of themedium with the scattering and absorption distributions being a multiplecomponent system, the total absorbance is given by the sum ofabsorbances of respective components.

[0019] The Inventor and others derived equations expressing therelations between various optical responses and optical characteristicsof scattering media from the MBL represented by above Eq (1) andclarified that it was possible to describe the various responses ofscattering media by an attenuation term dependent upon absorption and aterm independent of absorption, separated from each other. This made itfeasible to quantify the absolute value of the concentration of theabsorptive constituent in the scattering medium from the relationbetween the absorption coefficient and the extinction coefficientspecific to the absorptive constituent, based on spectral measurement ofthe attenuation term through use of light of multiple wavelengths.

[0020] The measuring methods based on the MBL as described above havethe major feature of being theoretically unsusceptible to the shape ofthe medium, the boundary conditions, and the scattering, and areapplicable to anisotropic scattering media and small media. However, theabsorption coefficient and the absorber concentration measured hereincan be correctly measured for scattering media with uniform absorption,but measurement for heterogeneous media with nonuniform absorptiondistribution provides a mean absorption coefficient and a meanconcentration (weighted mean for photon path distribution) ofphoton-passing portions. Accordingly, there arose the need fordetermining the photon path distribution, in order to realizehigh-accuracy optical CT.

[0021] Then the Inventor and others further expanded the methods ofanalyzing the behavior of photons in the scattering media on the basisof the MBL and have developed various measuring methods applicable evento the heterogeneous scattering media (heterogeneous system) withnonuniform distributions of scattering and absorption and capable ofquantitatively measuring a concentration distribution of an absorptiveconstituent in a scattering medium without being affected by the shapeof the medium, the scattering characteristics, and so on.

[0022] Specific measuring methods derived based on the MBL heretofore bythe Inventor and others include, for example, the time-resolvedspectroscopy (TRS) making use of the impulse response, thetime-integrated spectroscopy (TIS) making use of the time integral ofthe impulse response, the time-resolved gating integral spectroscopy(TGS) making use of the time-gating integral of the impulse response,the phase modulation spectroscopy (PMS) in the frequency domain, thereconstruction method of optical CT image based on the average valuemethod (AVM), and so on. These measuring methods are described inJapanese Patent Application No. 10-144300, aforementioned Document (7),and aforementioned Documents (15) to (25).

[0023] The materialization of the new measuring methods as describedabove drastically improved the measurement accuracy, but thereconstruction methods of optical CT image still used the conventionalweight function, i.e., the mean path length or the phase delayequivalent thereto. Specifically, a small change is given to anabsorption coefficient of voxel i in an imaginary medium with a uniformabsorption coefficient μ₈₄ described by the finite grid model, itsoptical output is calculated by making use of the Monte Carlocalculation, the numerical calculation of the Path Integral and thetransport equation, the numerical calculation of the photon diffusionequation, or the like, and the weight function (or the contributionfunction) Wi of voxel i is determined corresponding to the differencebetween optical outputs before and after the change, i.e., the deviationof light intensity, optical path length, attenuation, or the like withrespect to the imaginary medium having the uniform absorptioncoefficient μ_(ν).

DISCLOSURE OF THE INVENTION

[0024] However, the above-stated conventional measuring methods ofdetermining the weight function (or the contribution function)corresponding to the various measuring methods by changing theabsorption coefficient of voxel i in the imaginary medium described bythe finite grid model, required that the forward problem calculationsuch as the photon diffusion equation or the like should be executed thenumber of times equal to the total number of voxels, and thus had theproblem of extremely long calculation time and, in turn, the problem ofvery long measurement time. In this case, there was also the problem ofinsufficient measurement accuracy in the measuring methods using theconventional weight function in a narrow sense (the mean path length orthe phase delay equivalent thereto). For improvement in the measurementaccuracy, the calculation time became much longer, so as to result inthe problem of longer measurement time.

[0025] The present invention has been accomplished in view of theconventional problems as described above, and an object of the inventionis to provide a method that permits precise and quick calculation of adistribution of photon paths inside a scattering absorber of aheterogeneous system like living tissue. The photon path distributioncalculated in this way permits us to quickly calculate the weightfunction Wi corresponding to each of the various specific measuringmethods for quantitatively measuring the concentration distribution ofthe absorptive constituent inside the scattering medium.

[0026] In order to achieve the above object, the Inventor has conductedelaborate research on analytic methods and others based on the MBL, andfound that a time-resolved path length li in voxel i of a measuredmedium divided into N (1≦N) voxels was equal to a time-resolved pathlength li in voxel i (i=1, 2, . . . , N; 1≦N) of an imaginary mediumassumed to have the same scattering characteristics and boundaryconditions as the measured medium but absorb no light, that thistime-resolved path length li in voxel i of the imaginary medium could bedescribed by a photon existence probability or the like in voxel i (i=1,2, . . . , N; 1≦N) of the imaginary medium, and that the time-resolvedpath length li in voxel i of this imaginary medium could be directly andquickly calculated by making use of the result of calculation of photonmigration in the imaginary medium by the Monte Carlo calculation, thenumerical calculation of the path integral, the transport equation, orthe photon diffusion equation, or the like, thereby accomplishing thepresent invention.

[0027] The conventional Monte Carlo calculation is based on thefundamental theory of counting the number of occurrences of scatteringevents in voxel i in the finite grid model, and this concept is alsofollowed and is applied to the conventional weight function. On thisoccasion, if the scattering characteristics in voxel i are uniform, thepath length can be estimated from the relation between the number ofscattering events and the mean free pathlength. However, if thescattering characteristics are nonuniform in voxel i, the relationbetween the number of scattering events and the mean free pathlength isdependent upon a scattering distribution and becomes complicated, and itis thus infeasible to accurately determine the time-resolved path lengthbeing the objective of the present invention. Accordingly, in order tocircumvent the conventional problem as described, the present inventionemploys a new technique of calculating the time-resolved path length liin voxel i from the photon existence probability or the like in voxel iof the medium.

[0028] Specifically, a method of calculating a distribution of photonpaths inside a scattering medium according to the present invention is amethod comprising:

[0029] a first step of dividing a measured medium being the scatteringmedium, into N (1≦N) voxels each having a predetermined size and assumedto have a uniform absorption distribution in a measurement wavelengthregion;

[0030] a second step of hypothesizing an imaginary medium that can beassumed to have the same shape, boundary conditions, and scatteringdistribution as the measured medium, have a uniform refractive indexdistribution, and be nonabsorbing, for the measured medium divided intothe N voxels;

[0031] a third step of setting a light injection position for injectionof impulse light, and a light detection position for detection of animpulse response s(t) at a time t of the impulse light having beeninjected at the light injection position and having propagated in theimaginary medium, on a surface of the imaginary medium;

[0032] a fourth step of calculating the impulse response s(t) detectedat the light detection position after injection of the impulse light atthe light injection position into the imaginary medium;

[0033] a fifth step of calculating a photon existence probabilitydensity Ui(t′,t) that an ensemble of photons constituting the impulseresponse s(t) detected at the time t at the light detection positionafter injection of the impulse light at the light injection positioninto the imaginary medium has existed in an arbitrary voxel i (i=1, 2, .. . , N; 1≦N) in the imaginary medium at a time t′ (0≦t′≦t) before thedetection at the light detection position;

[0034] a sixth step of calculating a photon existence cumulativeprobability Ui(t) that the ensemble of photons constituting the impulseresponse s(t) detected at the time t at the light detection position hasexisted in said voxel i, using the photon existence probability densityUi(t′,t); and

[0035] an eighth step of calculating a time-resolved path length li inthe voxel i, using the photon existence cumulative probability Ui(t) andthe impulse response s(t).

[0036] Conventionally, a small change was given to the absorptioncoefficient in voxel i of the imaginary medium described by the finitegrid model and the various weight functions Wi were calculated from acalculated value of a ratio of intensities of output light before andafter the change of the absorption coefficient, a difference betweenattenuation amounts, or the like. This method required that the forwardproblem calculation of the photon diffusion equation or the like shouldbe executed the number of times equal to the total number N of voxels,and thus required the long calculation time, so as to result in theserious problem of very long measurement time, whereas the above methodsolved this problem. Moreover, since the time-resolved path length li iscalculated for the nonabsorbing medium in the above method, it alsosolved the problem that a long time was consumed for the forward problemcalculation of the photon diffusion equation or the like for theabsorbing medium in the conventional methods.

[0037] In the method of the present invention, parameters of all thevoxels are calculated in parallel for the set consisting of a pair ofthe light injection position and detection position, so that thecalculation time can be reduced more. Namely, the method of calculatingthe distribution of photon paths inside the scattering medium accordingto the present invention is the method first having permitted direct,precise, and quick calculation of the time-resolved path length li.

[0038] For specifically determining the time-resolved path length li invoxel i, as described hereinafter, thought is given to the photonexistence probability density Ui(t′,t) that the ensemble of photons withthe path length of l constituting the impulse response s(t) detected atthe time t at the light detection position after injection of theimpulse light at the light injection position into the finite gridmodel, has existed in the voxel i at the time t′ (0≦t′≦t) before thedetection at the light detection position.

[0039] The photon existence probability density Ui(t′,t) in the voxel iat the time t′, which is described by introducing the new time variablet′, can be expressed as follows by applying the reciprocity theorem.

U _(i)(t′,t)=U _(pi)(t′)U _(qi)(t−t′)   (2.2)

[0040] [In the above equation, Upi(t′) indicates a photon existenceprobability density that photons injected at the light injectionposition p at the time t′=0 exist in the voxel i at the time t′=t′, and

[0041] Uqi(t−t′) indicates a photon existence probability density thatphotons injected at the light detection position q at the time t−t′=0exist in the voxel i at the time t−t′.]

[0042] The above Ui(t′,t) can be determined by calculating two forwardproblems as indicated by Eq (2.2). At this time, Upi(t′) can becalculated by injection of the impulse light at the light injectionposition p of the finite grid model. Uqi(t−t′) can also be calculated byinjection of the impulse light at the light detection position q of thefinite grid model. A specific method of calculating the two forwardproblems represented by Eq (2.2) will be detailed later.

[0043] Since the present invention permits quick and direct calculationof the time-resolved path length li in the voxel i, it is furtherfeasible to quickly calculate the weighted mean Li of the time-resolvedpath length li with respect to time, using it.

[0044] Since there was no specific calculation method of thetime-resolved path length li known heretofore, no idea was proposedabout such an approach that the time-resolved path length li was firstdetermined and then the weighted mean Li of the time-resolved pathlength li and the weight function Wi were calculated based thereon.Actually, no report has been made heretofore on an example of themeasuring method of direct calculation and utilization of thetime-resolved path length li.

[0045] Since the present invention permits the quick and directcalculation of the time-resolved path length li in the voxel i, it isfurther feasible to quickly calculate the weight functions Wicorresponding to a variety of predetermined specific measuring methodsfor quantitatively measuring the concentration distribution of theabsorptive constituent inside the measured medium, using it.

[0046] Namely, the Inventor and others have already clarified therelations of the time-resolved path length li in the voxel i in themeasured medium with measured values by the various measuring methodsand the corresponding weight functions in the specific measuring methodsof TRS, TIS, TGS, PMS, AVM, etc. based on the MBL, in the case of use ofthe measured medium and the finite grid model resulting from thedivision of the corresponding imaginary medium into N voxels.

[0047] For example, the weight function for the attenuation amount isequivalent to the time-resolved path length li itself in TRS, and in TISit is described using the weighted mean Li (also referred to as the meanpath length Li) of the time-resolved path length li with respect totime, the variance of the time-resolved path length li, or the like.More specifically, the TIS and AVM making use of the attenuation amountand the mean path length employ the weight functions Wi described usingthe relations of the time-resolved path length li in the voxel i in themeasured medium with the mean path length Li being the weighted mean,the variance of the time-resolved path length li, etc., and it is knownthat these can be expressed by Eqs (7), (9), and (12) as describedhereinafter.

[0048] Here “time-resolved photon path distribution” in the method ofcalculating the photon path distribution inside the scattering mediumaccording to the present invention refers, as shown in FIG. 1, to aphoton path distribution of an ensemble of photons with the path lengthof l constituting an impulse response B measured at the light detectionposition q on the finite grid model 10 in the case where thethree-dimensional measured medium is assumed as the finite grid model 10divided into N (1≦N) voxels 20 (volume elements) and where impulse lightA is injected at the light injection position p on the finite grid model10. Specifically, this time-resolved photon path distribution isrepresented by an N-dimensional vector [li] having elements oftime-resolved path lengths li of the voxels i (i=1, 2, . . . , N; 1≦N).From the above definition, this time-resolved photon path distributionis a function of path length l, i.e., a function of time t. Theattenuation amount of output light signal in the various specificmeasuring methods is expressed by a function of the product of thetime-resolved path length li and the absorption coefficient of the voxeli. The weighted mean Li of the time-resolved path length li for time iscalled a mean path length in voxel i, and an N-dimensional vector [Li]having elements thereof of weighted means is called a mean path lengthdistribution. In the present specification, each “voxel 20” in FIG. 1and FIG. 2 will be denoted by “voxel i” instead of “voxel 20” forconvenience' sake of description.

[0049] Further, the weight function is defined as a contribution ofabsorption information to a physical quantity of interest. For example,the weight function for the time-resolved attenuation amount obtained inTRS is independent of absorption and absorption distribution, there isthe relation: [time-resolved attenuation amount]=[weightfunction]·[absorption coefficient difference], and this weight functionis equivalent to the time-resolved path length li obtained above. Theweight function for the attenuation amount obtained by the time integralof the impulse response (TRS) is defined by [TRS attenuationamount]=[weight function]·[absorption coefficient difference], and theweight function for the mean path length difference is defined by [meanpath length difference]=[weight function]·[absorption coefficientdifference]. In the case of TRS, these weight functions are functions ofabsorption coefficient.

[0050] By introducing these weight functions, the matrix equation forreconstruction of image in optical CT can be described as the product ofmatrices, i.e., as [weight function]·[absorption coefficientdifference], and the [absorption coefficient difference] or absorptiondistribution can be calculated from this equation.

[0051] The scattering characteristics of the foregoing measured mediumand imaginary medium can be heterogeneous, and the specific calculationof the forward problems can be performed utilizing the Monte Carlocalculation, the numerical calculation of the Path Integral or thetransport equation, the numerical calculation of the photon diffusionequation, and so on.

[0052] There are no specific restrictions on the “predeterminedmeasuring methods” to which the time-resolved path length li determinedby the method of calculating the photon path distribution inside thescattering medium according to the present invention is applied, and thepredetermined measuring methods can be, for example, the interiormeasuring methods for the scattering medium such as the time-resolvedspectroscopy (TRS), the time integral spectroscopy (TIS), thetime-resolved gating integral spectroscopy (TGS), the phase modulationspectroscopy in the frequency domain (PMS), the reconstruction method ofoptical CT image based on the average value method (AVM), and so onderived based on the MBL by the Inventor and others, and can also be themeasuring methods with no base on the MBL, as described in foregoingDocuments (1) to (12).

[0053] Another method of calculating a distribution of photon pathsinside a scattering medium according to the present invention is amethod comprising:

[0054] a first step of dividing a measured medium being the scatteringmedium, into N (1≦N) voxels each having a predetermined size and assumedto have a uniform absorption distribution in a measurement wavelengthregion;

[0055] a second step of hypothesizing an imaginary medium that can beassumed to have the same shape, boundary conditions, and scatteringdistribution as the measured medium, have a uniform refractive indexdistribution, and be nonabsorbing, for the measured medium divided intothe N voxels;

[0056] a third step of setting a light injection position for injectionof impulse light, and a light detection position for detection of animpulse response s(t) at a time t of the impulse light having beeninjected at the light injection position and having propagated in theimaginary medium, on a surface of the imaginary medium;

[0057] a fifth step of calculating a photon existence probabilitydensity Ui(t′,t) that an ensemble of photons constituting the impulseresponse s(t) detected at the time t at the light detection positionafter injection of the impulse light at the light injection positioninto the imaginary medium has existed in an arbitrary voxel i (i=1, 2, .. . , N; 1≦N) in the imaginary medium at a time t′ (0≦t′≦t) before thedetection at the light detection position;

[0058] a sixth step of calculating a photon existence cumulativeprobability Ui(t) that the ensemble of photons constituting the impulseresponse s(t) detected at the time t at the light detection position hasexisted in said voxel i, using the photon existence probability densityUi(t′,t);

[0059] a seventh step of calculating a sum U(t) of said photon existencecumulative probabilities Ui(t) for all the voxels; and

[0060] an eleventh step of calculating a time-resolved path length li inthe voxel i, using the photon existence cumulative probability Ui(t) andthe sum U(t) of the photon existence cumulative probabilities Ui(t) forall the voxels.

[0061] In the above case, the time-resolved path length li in the voxeli can be quickly and directly calculated, and thus the weighted mean Liof the time-resolved path length li with respect to time can be quicklycalculated using it. Since the time-resolved path length li in the voxeli can be quickly and directly calculated, it is feasible to quicklycalculate the weight functions Wi corresponding to the variouspredetermined specific measuring methods for quantitatively measuringthe concentration distribution of the absorptive constituent inside themeasured medium, using it.

BRIEF DESCRIPTION OF THE DRAWINGS

[0062]FIG. 1 is a schematic illustration showing a finite grid model foranalysis of photon migration inside a measured medium.

[0063]FIG. 2 is a schematic illustration showing paths of a photonensemble having passed the voxel i and a probability thereof among anensemble of photons (having the path length of l) constituting animpulse response B detected at the time t at the light detectionposition q in the case where the impulse light A is injected at thelight injection position p in the finite grid model shown in FIG. 1.

[0064]FIG. 3 is a graph showing an example of temporal waveforms of thephoton existence probability density Ui(t′,t) in voxel i at the time t′,the photon existence probability density Upi(t′) in voxel i at the timet′, and the photon existence probability density Uqi(t−t′) in voxel i atthe time t′ in the finite grid model shown in FIG. 2.

[0065]FIG. 4 is a graph schematically showing the relationship among thephoton existence cumulative probability Ui(t) in voxel i in the durationof 0≦t′≦t, the photon existence probability density Ui(t′,t) in voxel iat the time t′, and the photon existence probability density Upq(t′,t)in all the voxels at the time t′.

[0066]FIG. 5 is a flowchart showing an embodiment of the method ofcalculating the photon path distribution inside the scattering mediumaccording to the present invention.

[0067]FIG. 6 is a flowchart showing a second embodiment of the method ofcalculating the photon path distribution inside the scattering mediumaccording to the present invention.

BEST MODE FOR CARRYING OUT THE INVENTION

[0068] The preferred embodiments of the method of calculating the photonpath distribution inside the scattering medium according to the presentinvention will be described below in detail with reference to thedrawings. In the following description the same or equivalent portionswill be denoted by the same reference symbols and redundant descriptionwill be omitted. Although consideration using three-dimensionalcoordinates is necessary for photons traveling while being scattered inthe scattering medium, the following description will be sometimes givenusing two-dimensional coordinates for simplicity of description.

[0069] The principle of the present invention will be described first.

[0070] [Principle of the Invention]

[0071] As described previously, the analytic methods of analyzing thebehavior of light in the scattering medium, which have been developedbased on the MBL by the Inventor and others, can be applied to theheterogeneous scattering media (heterogeneous system) with nonuniformdistributions of scattering and absorption. It is then deduced from theanalytic methods based on the MBL that various optical responses (lightoutputs) from the heterogeneous scattering media can be described byseparating them into the attenuation term dependent upon absorption andthe term independent of absorption and that the attenuation term ofoptical response can be described by the photon path distribution andabsorption distribution in the measured medium. Accordingly, if thephoton path distribution in the measured medium can be known in advance,it will be feasible to implement optical CT or the like of quantifying aconcentration distribution of an absorptive constituent inside themeasured medium. Namely, important points for implementing the opticalCT are to know the photon path distribution inside the measured mediumand to develop a method of quickly calculating the photon pathdistribution.

[0072] The present invention is based on such starting points that theitems (i) to (v) below hold in the measured medium and in the finitegrid model of the corresponding imaginary medium.

[0073] Specifically, (i) supposing the refractive index distributioninside the measured medium is uniform, the relation of 1=ct always holdsfor the impulse response (time-resolved waveform) of the measuredmedium, where l is the path length, c the speed of photons in themeasured medium determined by the refractive index, and t the time. Thereason why the refractive index distribution in the measured medium issupposed to be uniform herein is that, where the heterogeneousscattering medium like living tissue is selected as a measured medium,the major ingredient of the living tissue is water and the refractiveindex distribution thereof can be assumed uniform. Therefore, thefollowing description will concern the measured medium with a uniformrefractive index distribution unless otherwise stated.

[0074] (ii) Since scattering and absorption are independent of eachother, the photon path distribution of an ensemble of photons (detectedas an impulse response) constituting an impulse response, i.e., thetime-resolved photon path distribution is uniquely determined once thelight injection position and light detection position are set for aheterogeneous medium with arbitrary distributions of scattering andabsorption.

[0075] (iii) The aforementioned photon path distribution is given by afunction of time, is independent of the number of injected photons(i.e., photon density) and of the absorption and absorption distributionof the measured medium, and is equivalent to a photon path distributionon the assumption that there occurs no absorption in the measuredmedium. When a three-dimensional measured medium is divided into N (1≦N)voxels (volume elements) each having a predetermined size, the photonpath distribution of the ensemble of photons having the path length of lcan be expressed by an N-dimensional vector [li] having elements of pathlengths li in the respective voxels i (i=1, 2, . . . , N; 1≦N).

[0076] (iv) This time-resolved photon path distribution [li] can becalculated by assuming an imaginary medium having the same scatteringcharacteristics (which can be heterogeneous) and boundary conditions asthe measured medium, having uniform absorption, and being nonabsorbing.

[0077] (v) The various weight functions Wi used for a variety of opticalCT image reconstruction algorithms can be described using thetime-resolved path length li in each voxel i defined above.

[0078] The Inventor proposed the optical CT for measuring theconcentration distribution of the absorptive constituent in theheterogeneous medium, using the weight function Wi described using thetime-resolved photon path distribution li and actually measured valuesfor the measured medium (Japanese Patent Application No. 10-144300 andDocument (19)).

[0079]FIG. 1 shows the finite grid model for analysis of photonmigration inside a heterogeneous scattering absorber as a measuredmedium.

[0080] Further, the finite grid model 10 shown in FIG. 1 also indicatesan imaginary medium (which will be denoted by the same reference symbolas the measured medium) assumed to have the same shape, scatteringdistribution, and boundary conditions as the measured medium, have auniform refractive index distribution, and be nonabsorbing, as well asthe measured medium. In this case, the three-dimensional measured mediumis divided into N (1≦N) voxels (volume elements) and the voxels aredenoted by respective numbers i (i=1, 2, . . . , N; 1≦N). Therefore, theimaginary medium assumed corresponding to the measured medium also has Nvoxels similarly.

[0081] It is noted herein that in the following description, when theabsorption coefficient of the medium represented by the finite gridmodel 10 in FIGS. 1 and 2 is zero, the “impulse response B” will besometimes described by designating it as “impulse response s(t),” i.e.,a response without absorption. When the absorption coefficient of themedium represented by the finite grid model 10 is not zero on the otherhand, the “impulse response B” will be sometimes described bydesignating it as “impulse response h(t),” i.e., a response withabsorption.

[0082] The size and shape of the voxel i in the finite grid model 10 ofthe imaginary medium shown in FIG. 1 are determined within the range andshape in which an absorption distribution in the voxel i of thecorresponding measured medium is assumed uniform or can be safelyassumed uniform. This is because the absorption distribution data ineach voxel i obtained by actual measurement is obtained as an average ofabsorption distributions in each voxel i. Once the size and the shape ofvoxel i are determined, the total number of all voxels is alsodetermined. Namely, the size and shape of each voxel i and the totalnumber of all voxels of the imaginary medium can be arbitrarilydetermined as long as the absorption can be assumed uniform in the voxeli, as described previously. The practical size and shape of voxel i aredetermined according to the measuring methods used, or according to theresolutions of CT image required by those measuring methods.

[0083] In the method of calculating the photon path distribution insidethe scattering medium according to the present invention, in the casewhere the impulse light A is injected at the light injection position pon the finite grid model 10 of the imaginary medium shown in FIG. 1, anensemble of photons having the path length of l is considered among theimpulse response B, i.e., h(t) measured at the light detection positionq of the finite grid model 10. Here the light injection position p andlight detection position q are arbitrary.

[0084] In the case where the impulse light (injected impulse) Aconsisting of a number of photons is injected into the finite grid model10 shown in FIG. 1, the impulse response B, or h(t) contains photonshaving various path lengths, but at a given observation time t the pathlength l is uniquely determined as 1=ct. There exist a lot of photonshaving the path length of l. Then let us define that the photon pathdistribution of the ensemble of photons having the path length of l,i.e., path lengths li in the respective voxels i (i=1, 2, . . . , N;1≦N) are represented by a probability distribution of photon paths thatthe photons can pass. In conjunction therewith, the photon pathdistribution in the voxels i in the finite grid model 10 is expressed byan N-dimensional vector [li] having elements of time-resolved pathlengths li in the voxels i. Since there exist a lot of photons havingthe path length of l, as described previously, the time-resolved pathlength li in each voxel i is an average (ensemble mean) of path lengthsof these individual photons in the voxel i. In the system as described,once the light injection position p and the light detection position qare determined, the time-resolved photon path distribution [li] of theensemble of photons constituting the impulse response B or h(t) isuniquely determined.

[0085] Since the time-resolved photon path distribution [li] of theensemble of photons constituting the impulse response h(t) of themeasured medium described above is independent of absorption of themeasured medium, it is equivalent to the time-resolved photon pathdistribution [li] of the ensemble of photons constituting the impulseresponse s(t) of the imaginary medium assumed to be nonabsorbing.Namely, by determining the time-resolved photon path distribution [li]of the ensemble of photons constituting the impulse response s(t) of theimaginary medium, it is feasible to determine the time-resolved photonpath distribution [li] of the ensemble of photons constituting theimpulse response h(t) of the measured medium.

[0086] 1. General Description of Analysis Model and Analysis

[0087] First described is the relation of the attenuation amount, pathlength, etc. in the measured medium being a heterogeneous scatteringmedium. Since the absorption and scattering phenomena are independent ofeach other, we can assume herein that the imaginary medium 10 shown inFIG. 1 has the same absorption distribution as the measured medium.However, the other conditions are the same as in the case of the finitegrid model 10 of the measured medium.

[0088] Then the relation between the path length l and the path lengthdistribution li is first given by the following equation from theaforementioned knowledges (i) to (v) and the MBL. $\begin{matrix}{{\sum\limits_{i = 1}^{N}l_{i}} = l} & (2)\end{matrix}$

[0089] [In the equation, i indicates each of numbers (i=1, 2, . . . , N;1≦N) of the N voxels in the measured medium, and li the time-resolvedpath length in voxel i.]

[0090] When the measured medium is a heterogeneous scattering medium,the impulse response h(t) is given by the following equationcorresponding to Eq (1). $\begin{matrix}{{h(t)} = {{s(t)}{\exp \lbrack {- {\sum\limits_{i = 1}^{N}( {\mu_{ai}l_{i}} )}} \rbrack}}} & (3) \\{\frac{{\partial\ln}\quad {h(t)}}{\partial\mu_{ai}} = {- l_{i}}} & (4) \\\begin{matrix}{{\ln \quad {h(t)}} = {{\ln \quad {s(t)}} - {\sum\limits_{i = 1}^{N}{\int_{0}^{\mu_{ai}}{l_{i}{\mu_{a}}}}}}} \\{= {{\ln \quad {s(t)}} - {\sum\limits_{i = 1}^{N}( {\mu_{ai}l_{i}} )}}}\end{matrix} & (5)\end{matrix}$

[0091] [In the equations, μ_(ai) indicates the absorption coefficient ofthe voxel i, s(t) the impulse response of the imaginary medium detectedat the detection position q of the imaginary medium on the assumption ofthe measured medium being nonabsorbing, and h(t) the impulse response ofthe measured medium detected at the detection position q. The last term(on the right side) in Eq (5) is derived from the fact that thetime-resolved path length li is independent of absorption and absorptiondistribution.]

[0092] Eq (4) and Eq (5) represent the differential form and theintegral form, respectively, of the impulse response h(t) represented byEq (3). It is seen from Eq (5) that the impulse response h(t) of theheterogeneous scattering medium can be described by separating it intothe term independent of absorption, lns(t), and the attenuation termdependent upon absorption, which is represented by Eq (16) below.$\begin{matrix}{\sum\limits_{i = 1}^{N}( {\mu_{ai}l_{i}} )} & (16)\end{matrix}$

[0093] Further, it is seen from Eq (5) that the attenuation amount ofthe optical response, ln[h(t)/s(t)], can be described by thetime-resolved photon path distribution [li] and the absorptiondistribution [μ_(ai)] in the measured medium. Attention has to be paidherein to the fact that the time-resolved photon path distribution inthe measured medium is equivalent to the time-resolved photon pathdistribution in the imaginary medium, i.e., in the nonabsorbing medium.In the time-resolved spectroscopy (TRS), therefore, the weight functionfor the attenuation amount is equal to the time-resolved path lengthdistribution [li], and it becomes feasible to implement reconstructionof optical CT image making use of this relation. Such time-resolvedoptical CT is substantiated by the easiest image reconstructionalgorithm.

[0094] More general optical CT includes the TIS making use of the timeintegral of the above impulse response h(t), the TGS making use of theintegral of the impulse response h(t) in a predetermined time range, thePMS making use of sinusoidal modulated light, and the AVM of providingthe imaginary medium with uniform absorption as a reference andmeasuring deviations from the reference, in those various methods. Forall of these methods, however, we can derive the fact that theattenuation amount or the deviation of attenuation amount can bedescribed by separating it into the attenuation term dependent uponabsorption and the term independent of absorption and the fact that theattenuation amount or the attenuation amount deviation of the variousoptical responses can be described by the time-resolved photon pathdistribution and the absorption distribution in the measured medium, anddetermine the weight function for the attenuation amount or theattenuation amount deviation from this relation. As describedhereinafter, it is also derived that the weight function for theattenuation amount or the attenuation amount deviation in the TIS andPMS can be described by the mean path length or the phase delay, andtheir higher-order partial differential coefficients.

[0095] For example, the optical response in the TIS is equivalent to thetime integral of the aforementioned impulse response h(t) and can beexpressed as in the following equation. $\begin{matrix}{I = {{\int_{0}^{\infty}{{h(t)}{t}}} = {\int_{o}^{\infty}{{s(t)}{\exp \lbrack {- {\sum\limits_{i = 1}^{N}( {\mu_{ai}l_{i}} )}} \rbrack}{t}}}}} & (6)\end{matrix}$

[0096] [In the equation, I indicates the time integral of the impulseresponse h(t) of the heterogeneous system.]

[0097] The differential form and integral form of Eq (6) can beexpressed as in the following equations, similar to Eq (4) and Eq (5)corresponding to Eq (3). $\begin{matrix}\begin{matrix}{\frac{{\partial\ln}\quad I}{\partial\mu_{ai}} = {- \frac{\int_{0}^{\infty}{l_{i}{s(t)}{\exp \lbrack {- {\sum\limits_{i = 1}^{N}( {\mu_{ai}l_{i}} )}} \rbrack}{t}}}{\int_{0}^{\infty}{{s(t)}{\exp \lbrack {- {\sum\limits_{i = 1}^{N}( {\mu_{ai}l_{i}} )}} \rbrack}{t}}}}} \\{= {{- {\langle l_{i}\rangle}} = {- L_{i}}}}\end{matrix} & (7) \\{{\ln \quad I} = {{\int_{0}^{\infty}{{s(t)}{t}}} - {\sum\limits_{i = 1}^{N}{\int_{0}^{\mu_{ai}}{{L_{i}( \mu_{a} )}{\mu_{a}}}}}}} & (8)\end{matrix}$

[0098] [In Eqs (7) and (8), Li indicates the mean path length Li invoxel i of an ensemble of photons constituting the time integral Idetected at the detection position, i.e., the weighted mean <li> of theaforementioned time-resolved path length li.]

[0099] This mean path length (weighted mean of the time-resolved pathlength li) Li is dependent upon the scattering characteristics, boundaryconditions, and absorption and absorption distribution of the scatteringmedium. On this occasion, attention has to be paid to the fact that theintegral in the second term on the right side in Eq (8) cannot beexpressed by the product, different from Eq (5). Accordingly, the weightfunction for the attenuation amount is not the simple mean path lengthLi, but is dependent upon absorption, which makes optical CTreconstruction algorithms complex.

[0100] In the AVM measurement of providing the imaginary medium withuniform absorption as a reference and measuring the measured medium asdeviations from the reference, the attenuation amount difference(attenuation amount deviation) ΔBi between the measured medium and theimaginary medium in voxel i is expressed as in the equation below, usingthe absorption coefficient difference (absorption coefficient deviation)Δμ_(ai) and the weight function Wi.

ΔB _(i))(μ_(ai))=Δμ_(ai) W _(i)   (9)

Δμ_(ai)=μ_(ai)−μ_(ν)  (10)

ΔB _(i)(μ_(ai))=B _(i)(μ_(ai))−B _(i)(μ₈₄)   (11)

[0101] $\begin{matrix}{W_{i} = {{L_{i}( \mu_{v} )} + {\frac{1}{2!}{\Delta\mu}_{ai}{L_{i}^{\prime}( \mu_{v} )}} + \ldots}} & (12)\end{matrix}$

B _(i)(μ_(ai))=∫₀ ^(μ) _(ai) L _(i)(μα)dμ _(a)   (13)

[0102] $\begin{matrix}{B = {\sum\limits_{i = 1}^{N}B_{i}}} & (14) \\\begin{matrix}{L = {{\langle l\rangle} = {\sum\limits_{i = 1}^{N}L_{i}}}} \\{= \frac{\int_{0}^{\infty}{{{ls}(t)}{\exp \lbrack {- {\sum\limits_{i = 1}^{N}( {\mu_{ai}l_{i}} )}} \rbrack}{t}}}{\int_{0}^{\infty}{{s(t)}{\exp \lbrack {- {\sum\limits_{i = 1}^{N}( {\mu_{ai}l_{i}} )}} \rbrack}{t}}}}\end{matrix} & (15)\end{matrix}$

[0103] [In Eqs (9) to (15), C indicates the fixed absorption coefficientof the imaginary medium,

[0104] Δμ_(ai) the absorption coefficient deviation of μ_(al) fromμ_(ν),

[0105] Bi (μ_(ai)) the attenuation amount in each voxel i about the timeintegral I of the impulse response h(t) of the measured medium,

[0106] Bi (μ_(ν)) the attenuation amount about the time integral I ofthe impulse response of the imaginary medium with the absorptioncoefficient being μ_(ν) (provided that this I is equal to a valueresulting when μ_(al)=μ_(ν) in Eq (6)),

[0107] ΔBi (μ_(al)) the attenuation amount deviation of Bi (μ_(al)) fromBi (μ_(ν)),

[0108] Li(μ_(ν)) the mean path length in voxel i of the imaginary mediumwith the absorption coefficient being μ_(ν) (the mean path length aboutthe time integral I, this Li (μ_(ν)) being equal to a value resultingwhen μ_(ai)=μ_(ν) in Eq (7)), which is equal to the weighted mean of thetime-resolved path length li in voxel i,

[0109] Wi the weight function for the attenuation amount deviationΔBi(μ_(ai)), which is represented by a function of the absorptioncoefficient deviation Δμ_(al) and first-order and higher-order partialdifferential coefficients (L, L′, L″, . . . ).]

[0110] Here the first-order and higher-order partial differentialcoefficients can be directly calculated from the time-resolved waveformdetermined for the imaginary medium.

[0111] The absorption coefficient distribution can also bequantitatively measured by using other various quantities derivedaccording to the various optical CT in addition to the above-statedmeasuring methods; for example, it can be measured by using acombination of the difference of mean path length with the weightfunction corresponding to the dispersion of path length, because arelation similar to above Eq (9) also holds between the difference ofmean path length and the weight function corresponding to the pathlength dispersion. Specifically, this measurement is performed for lightof different wavelengths and a concentration distribution of anabsorptive constituent (absolute value) is quantitatively determinedfrom the resultant absorption coefficients and the extinctioncoefficients specific to the substance.

[0112] 2. Method of Calculating Time-Resolved Path Length

[0113] Next, let us derive an equation for calculating the time-resolvedpath lengths li which are the elements of the time-resolved photon pathdistribution [li], for the imaginary medium described by the finite gridmodel 10.

[0114] Since it is feasible to calculate the various weight functionscorresponding to the measuring methods used in optical CT of the TIS(time-resolved integral spectroscopy), TGS (time-resolved gatingintegral spectroscopy), and PMS (phase modulation spectroscopy), fromthe time-resolved path length li, specific examples of these will alsobe described.

[0115] The voxel i of the measured medium has the absorption coefficientμ_(al), while the voxel i of the finite grid model 10 of the imaginarymedium has the absorption coefficient μ_(ai)=0. In this case, theimpulse response h(t) of the finite grid model 10 of the imaginarymedium shown in FIG. 1 can be expressed as follows by substitutingμ_(ai)=0 into Eq (3).

h(t)=s(t)   (1.1)

[0116] Namely, the impulse response s(t) of the imaginary mediumdetected at the light detection position q is equal to the impulseresponse h(t) of the measured medium detected at the light detectionposition q where the measured medium is assumed to be nonabsorbing.

[0117] At this time, the relation between the time-resolved path lengthli and the path length l can be expressed by the following equation.$\begin{matrix}{{\sum\limits_{i = 1}^{N}{l_{i}(t)}} = {l(t)}} & (2.1)\end{matrix}$

[0118] [In the equation, i indicates each of numbers i (i=1, 2, . . . ,N; 1≦N)) of the voxels constituting the finite grid model 10.]

[0119] The above relation between the time-resolved path length li andthe path length l also holds even in the case where a photon travels ina zigzag path to pass an identical voxel twice or more times. In theabove equation, the total path length l and the time-resolved pathlength li of each voxel i are represented by “l(t)” and “li(t),”respectively, in order to explicitly express that they are functions oftime t.

[0120] The following will detail the principle and the process ofcalculating the time-resolved path length li of voxel i, using such t asa parametric variable.

[0121]FIG. 2 is a schematic illustration showing the paths andprobability density of an ensemble of photons having passed the voxel iamong an ensemble of photons constituting the impulse response s(t) andhaving the path length of l (i.e., an ensemble of photons detected at apoint of time t), in the case where propagating light is detected at thelight detection position q after injection of impulse light A at thelight injection position p on the imaginary medium described by thefinite grid model 10 shown in FIG. 1, and the value Upiq(t) of such anensemble of photons detected at the point of time t at the lightdetection position q.

[0122] In the finite grid model 10, where propagating light is detectedat the light detection position q after injection of the pulse light Aat the light injection position p, let us consider a photon existenceprobability density Ui(t′,t) that an ensemble of photons constitutingthe impulse response s(t) and having the path length of l (i.e., anensemble of photons detected at the point of time t) exists in voxel iat a point of time t′ (0≦t′≦t) before the detection at the lightdetection position q. In the following description, this photonexistence probability density Ui(t′,t) will be referred to simply as“photon existence probability density Ui(t′,t) in voxel i at time t′.”

[0123] At this time, the photon existence probability density Ui(t′,t)in voxel i at time t′, described by introducing the aforementioned newtime variable t′, can be expressed like the following equation byapplying the reciprocity theorem.

U _(i)(t′,t)=U _(pi)(t′)U _(qi)(t−t′)   (2.2)

[0124] [In the equation, Upi(t′) indicates a photon existenceprobability density that photons injected at the light injectionposition p at time t′=0 exist in voxel i at time t′=t′,

[0125] Uqi(t−t′) indicates a photon existence probability density thatphotons injected at the light detection position q at time t−t′=0 existin voxel i at time t−t′.]

[0126] In the scattering media like living tissue, the reciprocitytheorem approximately holds for injection of collimated light andinjection of a pencil beam. For applying the reciprocity theorem moreprecisely, it is preferable to utilize a light injection method calledisotropic light injection. This is because the reciprocity theorem holdsmore accurately in the isotropic light injection method. This isotropiclight injection method is a light injection method reported by theInventor, which is disclosed, for example, in Japanese PatentApplication No. 5-301979 and Y. Tsuchiya, K. Ohta, and T. Urakami: Jpn.J. Appl. Phys. 34, Part 1, pp.2495-2501 (1995).

[0127] In the following description, Upi(t′) and Uqi(t−t′) in above Eq(2.2) will be referred to simply as “photon existence probabilitydensity Upi(t′) in voxel i at time t′” and as “photon existenceprobability density Upi(t−t′) in voxel i at time t−t′,” respectively.

[0128]FIG. 3 is a graph showing an example of the photon existenceprobability density Ui(t′,t) in voxel i at time t′ in the finite gridmodel 10, represented by above Eq (2.2), the photon existenceprobability density Upi(t′) in voxel i at time t′, and the photonexistence probability density Uqi(t−t′) in voxel i at time t−t′.

[0129] Above Ui(t′,t) can be determined by calculating the two forwardproblems indicated by Eq (2.2). At this time, Upi(t′) can be calculatedwith injection of the impulse light A at the light injection position pon the finite grid model 10. Uqi(t−t′) can also be calculated withinjection of the impulse light A at the light detection position q onthe finite grid model 10. For specific calculation of the above forwardproblems, it is possible to apply the Monte Carlo calculation, thenumerical calculation of Path Integral or the transport equation, thenumerical calculation of the photon diffusion equation, and so on. Onthis occasion, since the finite grid model 10 is nonabsorbing, thecalculation algorithm becomes simpler than the calculation of forwardproblems in the conventional absorbing case, so that the calculationtime also becomes shorter.

[0130] By calculating the aforementioned photon existence probabilitydensity Ui(t′,t) in voxel i at time t′, it is feasible to calculateaccording to the equation below, the probability density Ui(t) that anensemble of photons with the path length of l constituting the impulseresponse s(t) detected at the point of time t at the light detectionposition q after injection of the impulse light A at the light injectionposition p into the finite grid model 10 exists in voxel i during theperiod of 0≦t′≦t in which photons travel from the light injectionposition p to the light detection position q in the finite grid model10. This probability density Ui(t) is accumulation of probabilitiesincluding one that a photon traveling in a zigzag path in the finitegrid model 10 (in the imaginary medium or in the measured medium) hasexisted in an identical voxel at different times.

∫₀ ^(t) U _(i)(t′,t)dt′=U _(i)(t)   (2.3)

[0131] Ui(t) represented by above Eq (2.3) indicates a photon existencecumulative probability in voxel i of an ensemble of photons constitutingthe impulse response s(t) and having the path length of l, and in thefollowing description, it will be referred to simply as “photonexistence cumulative probability Ui(t) in voxel i during the period0≦t′≦t.”

[0132] Here the summation of the photon existence cumulative probabilityUi(t) represented by Eq (2.3), over all voxels in the finite grid model10 provides a photon existence cumulative probability U(t) in the finitegrid model 10 to yield the following equation. $\begin{matrix}{{\sum\limits_{i = 1}^{N}{\int_{0}^{t}{{U_{i}( {t^{\prime},t} )}{t^{\prime}}}}} = {{\sum\limits_{i = 1}^{N}{U_{i}(t)}} = {U(t)}}} & (2.4)\end{matrix}$

[0133] This photon existence cumulative probability U(t) indicates acumulative probability that an ensemble of photons with the path lengthof l constituting the impulse response s(t) detected at the point oftime t at the light detection position q after injection of the impulselight A at the light injection position p into the finite grid model 10exists in the finite grid model 10 (in all the voxels) during the period0≦t′≦t in which the photons travel in the finite grid model 10 from thelight injection position p to the light detection position q.

[0134] In the following description, U(t) represented by above Eq (2.4)will be referred to as “photon existence cumulative probability U(t) inall voxels during period 0≦t′≦t.”

[0135] On the other hand, against the aforementioned photon existenceprobability density Ui(t′,t) in voxel i at time t′, let us consider aprobability density Upq(t′,t) that an ensemble of photons with the pathlength of l constituting the impulse response s(t) detected at the pointof time t at the light detection position q after injection of theimpulse light A at the light injection position p into the finite gridmodel 10 exists in the finite grid model 10 (in the imaginary medium) atthe point of time t′ (0≦t′≦t) prior to the detection at the lightdetection position q. Namely, let us consider a “photon existenceprobability density Upq(t′,t) in all voxels at time t′,” against the“photon existence probability density Ui(t′,t) in voxel i at time t′.”In this case, since consideration is given to the photon existenceprobability density at time t′, the summation of photon existenceprobability densities over all voxels is equal to Upq(t′,t) and thisUpq(t′,t) can be expressed by the following equation using Ui(t′,t).$\begin{matrix}{{U_{pq}( {t^{\prime},t} )} = {\sum\limits_{i = 1}^{N}{U_{i}( {t^{\prime},t} )}}} & (2.5)\end{matrix}$

[0136] The definite integral of this Eq (2.5) over the range from thetime t′=0 to the time t′=t is naturally equal to the previously obtainedphoton existence cumulative probability U(t) in all the voxels duringthe period 0≦t′≦t, and can be expressed by the following equation.

U(t)=∫₀ ^(t) U _(pq)(t′,t)dt′  (2.6)

[0137] Since the nonabsorbing imaginary medium is assumed in the methodof calculating the photon path distribution inside the scattering mediumaccording to the present invention herein, the photon existenceprobability density Upq(t′,t) in all voxels at time t′ is constant fromthe time t′=0 to time t′=t. Then the cumulative probability U(t) can beexpressed like the following equation using Eq (2.6).

U(t)=tU _(pq)(t′,t)   (2.7)

[0138] In addition to the above fact, since the photon existenceprobability density Upq(t′,t) in all voxels at time t′ in thenonabsorbing imaginary medium is constant from the time t′=0 to the timet′=t, it is understood that this Upq(t′,t) is equal to the impulseresponse detected at the point of time t at the light detection positionq, i.e., s(t) defined in Eq (1.1). Namely, the relationship betweenUpq(t′,t) and s(t) is expressed like the following equation.

s(t)=U _(pq)(t′,t)   (2.8)

[0139] From the above, the photon existence cumulative probability U(t)in all voxels during period 0≦t′≦t in Eq (2.4) is eventually expressedby the following equation, by substituting the result of Eq (2.8) intoEq (2.7).

U(t)=ts(t)   (2.9)

[0140] Further, it is understood that the summation of aforementionedUpiq(t) described with FIG. 2, over all voxels is equal to the photonexistence cumulative probability U(t) in the finite grid model 10 andthe relationship indicated by the following equation holds.$\begin{matrix}{{\sum\limits_{i = 1}^{N}{U_{piq}(t)}} = {{\sum\limits_{i = 1}^{N}{U_{i}(t)}} = {U(t)}}} & (2.10)\end{matrix}$

[0141]FIG. 4 schematically shows the relationship among theaforementioned photon existence cumulative probability Ui(t) in voxel iduring period 0≦t′≦t, the photon existence probability density Ui(t′,t)in voxel i at time t′, and the photon existence probability densityUpq(t′,t) in all voxels at time t′. Here reference numeral 30 in FIG. 4designates a photon.

[0142] Let li(t′,t) be the time-resolved path length in voxel i at timet′ of an ensemble of photons with the path length of l among photonsconstituting the impulse response s(t) detected at the point of time tat the light detection position q after injection of the pulse light Aat the light injection position p into the finite grid model 10. Thenthe integral of this li(t′,t) from the time t′=0 to the time t′=t isequal to the time-resolved path length li(t) and can be expressed by thefollowing equation.

∫₀ ^(t) li(t′,t)dt′=l _(i)(t)   (2.11)

[0143] Incidentally, since the speed of photons in the finite grid model10 is constant, the path length is always proportional to the photonexistence probability at an arbitrary time or an arbitrary portion.

[0144] Specifically, when focus is placed on an ensemble of photons withthe path length of l among photons constituting the impulse response s(t) detected at the point of time t at the light detection position qafter injection of the impulse light A at the light injection position pinto the finite grid model 10, the time-resolved path length li(t′,t) invoxel i at an arbitrary time t′ is proportional to the photon existenceprobability density Ui(t′,t) in voxel i at the time t′. Then thetime-resolved path length li(t), which is the time integral of li(t′,t),is also proportional to the photon existence cumulative probabilityUi(t) which is the time integral of Ui(t′,t). Further, when focus isplaced on an ensemble of photons with the path length of l among photonsconstituting the impulse response s(t) detected at the point of time tat the light detection position q after injection of the impulse light Aat the light injection position p into the finite grid model 10, thepath length l is proportional to the photon existence cumulativeprobability U(t).

[0145] By letting A be a constant of the proportionality, the followingequation can be yielded from the relations of Eq (2.1), Eq (2.9), and Eq(2.11). $\begin{matrix}{A = {\frac{l_{i}( {t^{\prime},t} )}{U_{i}( {t^{\prime},t} )} = {\frac{l_{i}(t)}{U_{i}(t)} = {\frac{\sum\limits_{i = 1}^{N}{l_{i}(t)}}{\sum\limits_{i = 1}^{N}{U_{i}(t)}} = {\frac{l(t)}{U(t)} = {\frac{l(t)}{{ts}(t)} = \frac{c}{s(t)}}}}}}} & (2.12)\end{matrix}$

[0146] [In the equation, c is the speed of photons traveling in themedium.]

[0147] By performing equivalent deformation of the relation of Eq (2.12)and solving it with respect to the time-resolved path length li(t), weobtain the following equation. $\begin{matrix}{{l_{i}(t)} = {{\frac{U_{i}(t)}{U(t)}{l(t)}} = {{\frac{U_{i}(t)}{{ts}(t)}{l(t)}} = {\frac{c}{s(t)}{U_{i}(t)}}}}} & (2.13)\end{matrix}$

[0148] From this Eq (2.13), the elements expressing the time-resolvedphoton path distribution as the N-dimensional vector [li], i.e., thetime-resolved path length li for each voxel can be calculated bydetermining Ui(t′,t) by calculating the forward problems indicated by Eq(2.2) and further using Ui(t) calculated from this Ui(t′,t), U(t), and,l determined as 1=ct by given time t, or s(t).

[0149] The time-resolved path length li is directly used in thetime-resolved optical CT image reconstruction making use of thetime-resolved response. This is because the weight function intime-resolved optical CT is equal to the time-resolved path length li.As described hereinafter, the time-resolved path length li is the basequantity for describing the weight functions in the various measuringmethods.

[0150] 3. Specific Calculation Methods of Time-Resolved Path Length li

[0151]FIG. 5 is a flowchart showing an embodiment of the method ofcalculating the photon path distribution inside the scattering mediumaccording to the present invention. The following will describe theflowchart shown in FIG. 5.

[0152] The first step is to divide the measured medium as a scatteringmedium into N (1≦N) voxels each having a predetermined size and assumedto have a uniform absorption distribution in a measurement wavelengthregion (S1).

[0153] The second step is to hypothesize an imaginary medium that can beassumed to have the same shape, boundary conditions, and scatteringdistribution as the measured medium, have a uniform refractive indexdistribution, and be nonabsorbing, for the measured medium divided intothe N voxels (S2).

[0154] The third step is to set, on a surface of the imaginary medium, alight injection position p for injection of impulse light A, and a lightdetection position q for detection of an impulse response s(t) at a timet of the impulse light A having been injected at the light injectionposition p and having propagated in the imaginary medium (S3).

[0155] The fourth step is to calculate the impulse response s(t)detected at the light detection position after injection of the impulselight A at the light injection position p into the imaginary medium(S4).

[0156] The next step is to calculate the photon existence probabilitydensity Upi(t′) that, after injection of the impulse light A at thelight injection position p into the imaginary medium, photons injectedat the light injection position p at the time t′=0 exist in voxel i atthe time t′=t′. Subsequently calculated is the photon existenceprobability density Uqi(t−t′) that, after injection of the impulse lightA at the light detection position q into the imaginary medium, photonsinjected at the light injection position p at the time t−t′=0 exist invoxel i at the time t−t′. Then calculated according to Eq (2.2) is thephoton existence probability density Ui(t′,t) that, after injection ofthe impulse light A at the light injection position p into the finitegrid model 10, an ensemble of photons with the path length of lconstituting the impulse response s(t) detected at the point of time tat the light detection position exists in voxel i at the point of timet′ (0≦t′≦t) prior to the detection at the light detection position (S5).

[0157] In the above case, the photon existence probability densities inall the voxels constituting the finite grid model 10 can be calculatedin parallel or simultaneously for one light injection position. For thisreason, there is no need for the conventional time-consuming calculationof varying the absorption coefficient of each voxel and repeatedlysolving the forward problems by the number of times equal to the totalnumber N of voxels, which extremely decreases the calculation time.

[0158] The next step is to calculate the photon existence cumulativeprobability Ui(t) that the ensemble of photons constituting the impulseresponse s(t) detected at the point of time t at the light detectionposition q has existed in voxel i, from the photon existence probabilitydensity Ui(t′,t) according to Eq (2.3) (S6).

[0159] The last step is to calculate the time-resolved path length li invoxel i, using the photon existence cumulative probability Ui(t) and theimpulse response s(t) according to Eq (2.13) (S8).

[0160] As described in the first embodiment of the present inventionabove, the method of calculating the photon path distribution inside thescattering medium according to the present invention is configured tosolve the forward problems indicated by Eq (2.2) for the imaginarymedium assumed to have the same scattering characteristics and boundaryconditions as the measured medium, and be nonabsorbing, and thus thecalculation time thereof becomes shorter than in the conventionalabsorbing case of computation of forward problems. Further, it obviatesthe need for the conventional time-consuming problem calculation ofvarying the absorption coefficient of each voxel and repeatedly solvingthe forward problems by the number of times equal to the total number ofvoxels, and permits direct and quick calculation of the time-resolvedpath length li in voxel i.

[0161] When the above method of calculating the photon path distributioninside the scattering medium is applied to the optical CT using variouslight responses, the weight function Wi corresponding to thepredetermined measuring method for quantitatively measuring theconcentration distribution of the absorptive constituent inside themeasured medium is calculated using the time-resolved path length licalculated for each arbitrary voxel i of the imaginary medium in thelast stage of the flowchart shown in FIG. 5. In this case, there are amethod of directly calculating the weight function Wi corresponding tothe predetermined measuring method from the time-resolved path length licalculated for each arbitrary voxel i of the imaginary medium; and amethod of first determining the weighted mean Li for time of thetime-resolved path length li calculated for each arbitrary voxel i ofthe imaginary medium and then calculating the weight function Wicorresponding to the predetermined measuring method using the weightedmean Li.

[0162] As also described previously, there was no conventionally knownspecific calculation method of the time-resolved path length li and,therefore, there was no idea proposed on the method of first determiningthe time-resolved path length li and thereafter calculating the weightedmean Li of the time-resolved path length li and the weight function Wi.In contrast to it, since the calculation methods of the weight functionherein are based on the quick and direct calculation of thetime-resolved path length li of voxel i, they permit quick calculationof the weighted mean Li of the time-resolved path length li or theweight function Wi corresponding to each of the various predeterminedspecific measuring methods for quantitatively measuring theconcentration distribution of the absorptive constituent in the measuredmedium, using the time-resolved path length li. Specific calculationmethods of the weight functions corresponding to optical CT usingvarious light responses, using the time-resolved path length licalculated for each arbitrary voxel i of the imaginary medium will bedescribed later.

[0163]FIG. 6 is a flowchart showing a second embodiment of the method ofcalculating the photon path distribution inside the scattering mediumaccording to the present invention. The following will describe theflowchart shown in FIG. 6.

[0164] The first step is to divide the measured medium as a scatteringmedium into N (1≦N) voxels each having a predetermined size and assumedto have a uniform absorption distribution, in the measurement wavelengthregion (S1).

[0165] The second step is to hypothesize an imaginary medium that can beassumed to have the same shape, boundary conditions, and scatteringdistribution as the measured medium, have a uniform refractive indexdistribution, and be nonabsorbing, for the measured medium divided intothe N voxels (S2).

[0166] The third step is to set, on a surface of the imaginary medium, alight injection position p for injection of impulse light A, and a lightdetection position q for detection of the impulse response s(t) at thetime t of the impulse light A having been injected at the lightinjection position p and having propagated in the imaginary medium (S3).

[0167] The next step is to calculate the photon existence probabilitydensity Upi(t′) that, after injection of the impulse light A at thelight injection position p into the imaginary medium, photons injectedat the light injection position p at the time t′=0 exist in voxel i atthe time t′=t′. Subsequently calculated is the photon existenceprobability density Uqi(t−t′) that, after injection of the impulse lightA at the light detection position q into the imaginary medium, photonsinjected at the light injection position p at the time t−t′=0 exist invoxel i at the time t−t′. Then calculated according to Eq (2.2) is thephoton existence probability density Ui(t′,t) that, after injection ofthe impulse light A at the light injection position p into the finitegrid model 10, an ensemble of photons with the path length of lconstituting the impulse response s(t) detected at the point of time tat the light detection position exist in voxel i at the point of time t′(0≦t′≦t) prior to the detection at the light detection position (S5).

[0168] In the above case, the photon existence probability densities inall the voxels constituting the finite grid model 10 can also becalculated in parallel or simultaneously for one light injectionposition. For this reason, there is no need for the conventionaltime-consuming calculation of varying the absorption coefficient of eachvoxel and repeatedly solving the forward problems by the number of timesequal to the total number N of voxels, which extremely decreases thecalculation time.

[0169] The next step is to calculate the photon existence cumulativeprobability Ui(t) that the ensemble of photons constituting the impulseresponse s(t) detected at the point of time t at the light detectionposition q has existed in voxel i, from the photon existence probabilitydensity Ui(t′,t) according to Eq (2.3) (S6).

[0170] The next step is to calculate a sum U(t) of the photon existencecumulative probabilities Ui(t) for all the voxels according to Eq (2.4)(S7).

[0171] The last step is to calculate the time-resolved path length li invoxel i, using the photon existence cumulative probability Ui(t) and thesum U(t) over all the voxels according to Eq (2.13) (S11).

[0172] It is seen from the above description of the second embodiment ofthe present invention, as in the description of the first embodiment ofthe present invention described previously, that the method ofcalculating the photon path distribution inside the scattering mediumaccording to the present invention is configured to solve the forwardproblems indicated by Eq (2.2) for the imaginary medium assumed to havethe same scattering characteristics and boundary conditions as themeasured medium and be nonabsorbing and thus the calculation timethereof becomes shorter than in the conventional absorbing case ofcomputation of forward problems. It is further seen that it obviates theneed for the conventional time-consuming problem calculation of varyingthe absorption coefficient of each voxel and repeatedly solving theforward problems by the number of times equal to the total number ofvoxels and thus permits direct and quick calculation of thetime-resolved path length li in voxel i.

[0173] When the above method of calculating the photon path distributioninside the scattering medium is applied to the optical CT using variouslight responses, as in the first embodiment of the present inventiondescribed previously, the weight function Wi corresponding to thepredetermined measuring method for quantitatively measuring theconcentration distribution of the absorptive constituent inside themeasured medium is calculated using the time-resolved path length licalculated for each arbitrary voxel i of the imaginary medium in thelast stage of the flowchart shown in FIG. 6. In this case, there arealso the method of directly calculating the weight function Wicorresponding to the predetermined measuring method from thetime-resolved path length li calculated for each arbitrary voxel i ofthe imaginary medium; and the method of first determining the weightedmean Li for time of the time-resolved path length li calculated for eacharbitrary voxel i of the imaginary medium and then calculating theweight function Wi corresponding to the predetermined measuring methodusing the weighted mean Li.

[0174] Since the calculation methods of weight function are based on thequick and direct calculation of the time-resolved path length li ofvoxel i, they permit quick calculation of the weighted mean Li of thetime-resolved path length li or the weight function Wi corresponding toeach of the various predetermined specific measuring methods forquantitatively measuring the concentration distribution of theabsorptive constituent inside the measured medium, using thetime-resolved path length li.

[0175] As described above, the methods of calculating the photon pathdistribution inside the scattering medium according to the presentinvention can be applied to the optical CT using various lightresponses. In this case, the weight function Wi corresponding to thepredetermined measuring method for quantitatively measuring theconcentration distribution of the absorptive constituent inside themeasured medium is calculated using the time-resolved path length licalculated for each arbitrary voxel i of the imaginary medium, in thelast stage of the flowchart shown in FIG. 5 or in FIG. 6.

[0176] For example, specific examples include the time-resolvedspectroscopy (TIS), the time-resolved gating integral spectroscopy(TGS), the phase modulation spectroscopy in the frequency domain (PMS),the average value method (TIS), etc. derived based on the MBL by theInventor and others. In these cases, various light responses are definedby mathematical conversion of the aforementioned impulse response, andthe optical CT processes corresponding to the respective methods employtheir respective weight functions different from each other.

[0177] The following will describe examples of application of thetime-resolved path length li determined by the method of calculating thephoton path distribution inside the scattering medium according to thepresent invention to the optical CT using the light responses in thetime-resolved spectroscopy (TIS), the time-resolved gating integralspectroscopy (TGS), the phase modulation spectroscopy in the frequencydomain (PMS), and so on.

[0178] 4. Photon Path Distribution Inside Scattering Medium and WeightFunction in TIS

[0179] As described previously, the optical CT based on TIS using thetime integral I of the impulse response h(t) employs the weight functiondescribed using the mean path length Li of voxel i or the path lengthdispersion or the like for the impulse response h(t). This weightfunction Wi is defined as a weight function for a reference medium witha uniform absorption distribution.

[0180] First, let us consider the time integral I of the impulseresponse h(t) of a scattering medium with a heterogeneous absorptiondistribution. A model considered herein is one obtained by setting theabsorption coefficient of voxel i of the imaginary medium described bythe finite grid model 10 shown in FIG. 1, as μ_(al). Using these, thetime integral I of the impulse response h(t) is given by the followingequation. $\begin{matrix}{{\ln \quad I} = {{\ln \quad {\int_{0}^{\infty}{{s(t)}{t}}}} - {\sum\limits_{i = 1}^{N}{\int_{0}^{\mu_{a_{i}}}{{L_{i}( \mu_{a} )}{\mu_{a}}}}}}} & (4.1)\end{matrix}$

[0181] [In the equation, μ_(ai) indicates the absorption coefficient ofvoxel i,

[0182] s(t) indicates the impulse response h(t) detected at thedetection position on the assumption that the measured medium isnonabsorbing, and

[0183] Li indicates the mean path length in voxel i (weighted mean oftime-resolved path length li in voxel i) of an ensemble of photonsconstituting the time integral I detected at the detection position.]

[0184] In this equation, InI is described as separated into the termindependent of absorption and the term dependent upon the absorptioncoefficient μ_(ai) of each voxel (attenuation).

[0185] Here the mean path length Li of voxel i is a quantity dependentupon the scattering characteristics, boundary conditions, and absorptionand absorption distribution of the scattering medium, and is given bythe following equation. $\begin{matrix}{\frac{{\partial\ln}\quad I}{\partial\mu_{ai}} = {{\frac{1}{I}\frac{\partial I}{\partial\mu_{ai}}} = {\frac{\int_{0}^{\infty}{l_{i}{s(t)}{\exp \lbrack {- {\sum\limits_{i = 1}^{N}( {\mu_{ai}l_{i}} )}} \rbrack}{t}}}{\int_{0}^{\infty}{{s(t)}{\exp \lbrack {- {\sum\limits_{i = 1}^{N}( {\mu_{ai}l_{i}} )}} \rbrack}{t}}} = {- L_{i}}}}} & (4.2)\end{matrix}$

[0186] The weight function Wi used in the corresponding optical CT isdefined for the imaginary medium with the uniform absorption coefficientμ_(ν), and is given by the following equation. A model in this case isone obtained by setting the absorption coefficient of voxel i in themodel shown in FIG. 1, as μ_(ν). $\begin{matrix}{W_{i} = {{L_{i}( \mu_{v} )} + {\frac{1}{2!}{\Delta\mu}_{ai}{L_{i}^{\prime}( \mu_{v} )}} + \ldots}} & (4.3)\end{matrix}$

[0187] Δμ_(ai)=μ_(ai)−μ_(ν)  (4.4) $\begin{matrix}{{L( \mu_{v} )} = {\sum\limits_{i = 1}^{N}{L_{i}( \mu_{v} )}}} & (4.5) \\{{L_{i}( \mu_{v} )} = \frac{\int_{0}^{\infty}{l_{i}{s(t)}{\exp \lbrack {{- \mu_{v}}l} \rbrack}{t}}}{\int_{0}^{\infty}{{s(t)}{\exp \lbrack {{- \mu_{v}}l} \rbrack}{t}}}} & (4.6)\end{matrix}$

[0188] [In Eqs (4.3) to (4.6), μ_(ν) indicates the constant absorptioncoefficient given to the imaginary medium,

[0189] Δμ_(ai) the absorption coefficient deviation of μ_(ai) fromμ_(ν),

[0190] Li (μ_(ν)) the mean path length in voxel i (weighted mean oftime-resolved path length li in voxel i) of the imaginary medium withthe constant absorption coefficient μ_(ν),

[0191] L (μ_(ν)) the mean path length of the impulse response h(t) ofthe imaginary medium with the constant absorption coefficient μ_(ν) (theweighted mean of the path length l in the imaginary medium with theabsorption coefficient μ_(ν)), and

[0192] Wi the weight function for the attenuation amount deviation ΔBi(μ_(ai)) in aforementioned Eq (9), which is represented by a function ofthe absorption coefficient deviation Δμ_(ai) and the first-order andhigher-order partial differential coefficients (L, L′, L″, . . . ).]

[0193] Here the first-order and higher-order partial differentialcoefficients can be directly calculated from the time-resolved waveformobtained for the imaginary medium.

[0194] 5. Photon Path Distribution Inside Scattering Medium and WeightFunction in TGS

[0195] First, a time-resolved gating integral signal Ig of the impulseresponse h(t) of a scattering medium with a heterogeneous absorptiondistribution is given by the following equation. $\begin{matrix}{{\ln \quad I_{g}} = {{\ln \quad {\int_{t_{1}}^{t_{2}}{{s(t)}{t}}}} - {\sum\limits_{i = 1}^{N}{\int_{0}^{\mu_{ai}}{{L_{gi}( \mu_{a} )}{\mu_{a}}}}}}} & (5.1)\end{matrix}$

[0196] [In the equation, L_(gl) indicates a mean path length in voxel i(weighted mean of time-resolved path length li in a gate time in voxeli) of an ensemble of photons constituting the time gating integral Ig ofthe impulse response h(t).]

[0197] Here the time range of the gating integration is [t1,t2], and bysetting it as [0, ∞], the gating integral yields the time integralsignal I of the impulse response h(t) obtained in the previous section.

[0198] The weight function Wgi used in the corresponding optical CT isdefined for the imaginary medium with the uniform absorption coefficientμ_(ν), and is given by the following equation. $\begin{matrix}{W_{gi} = {{L_{qi}( \mu_{v} )} + {\frac{1}{2!}{\Delta\mu}_{ai}{L_{gi}^{\prime}( \mu_{v} )}} + \ldots}} & (5.2) \\{{L_{g}( \mu_{v} )} = {\sum\limits_{i = 1}^{N}{L_{gi}( \mu_{v} )}}} & (5.3) \\{{L_{gi}( \mu_{v} )} = \frac{\int_{t_{1}}^{t_{2}}{l_{i}{s(t)}{\exp \lbrack {{- \mu_{v}}l} \rbrack}{t}}}{\int_{t_{1}}^{t_{2}}{{s(t)}{\exp \lbrack {{- \mu_{v}}l} \rbrack}{t}}}} & (5.4)\end{matrix}$

[0199] [In Eqs (5.2) to (5.4), μ_(ν) indicates the constant absorptioncoefficient given to the imaginary medium,

[0200] Δμ_(al) indicates the absorption coefficient deviation of μ_(ai)from μ_(ν), as presented in aforementioned Eq (4.4),

[0201] Lgi (μ_(ν)) indicates the mean path length in voxel i (weightedmean of time-resolved path length li in gating time in voxel i) of anensemble of photons constituting the time gating integral Ig of theimpulse response h(t) of the imaginary medium,

[0202] Lg (μ_(ν)) indicates the mean path length (weighted mean) of anensemble of photons constituting the time gating integral Ig of theimpulse response h(t) of the imaginary medium given the constantabsorption coefficient μ_(ν), and

[0203] Wi indicates the weight function for the attenuation amountdeviation ΔBgi (μ_(ai)) in application of the attenuation amountdeviation ΔBi (μ_(ai)) in aforementioned Eq (9) to TGS, which isrepresented by a function of the absorption coefficient deviationΔμ_(ai) and the first-order and higher-order partial differentialcoefficients (Lg, Lg′, Lg″, . . . ).]

[0204] 6. Photon Path Distribution Inside Scattering Medium and WeightFunction in PMS

[0205] For a scattering medium with a heterogeneous absorptiondistribution, a response obtained in PMS is the Fourier transform of theimpulse response h(t) for the scattering medium with the heterogeneousabsorption distribution and is given by the following equation.$\begin{matrix}\begin{matrix}{{H(\omega)} = {\int_{0}^{\infty}{{s(t)}{\exp ( {- b} )}{\exp ( {{- {j\omega}}\quad t} )}{t}}}} \\{= {{R( {b,\omega} )} + {j\quad {X( {b,\omega} )}}}} \\{= {{A( {b,\omega} )}{\exp \lbrack {- {{j\varphi}( {b,\omega} )}} \rbrack}}}\end{matrix} & (6.1) \\{b = {{\sum\limits_{i = 1}^{N}( {\mu_{ai}l_{i}} )} = {\sum\limits_{i = 1}^{N}( {\mu_{ai}{ct}_{i}} )}}} & (6.2)\end{matrix}$

[0206] Here R, X, A, and φ are defined as follows. $\begin{matrix}{R = {\int_{0}^{\infty}{{s(t)}{\exp \lbrack {- {\sum\limits_{i = 1}^{N}( {\mu_{ai}{ct}_{i}} )}} \rbrack}\cos \quad \omega \quad t{t}}}} & \text{(6.3)} \\{X = {- {\int_{0}^{\infty}{{s(t)}{\exp \lbrack {- {\sum\limits_{i = 1}^{N}( {\mu_{ai}{ct}_{i}} )}} \rbrack}\sin \quad \omega \quad t{t}}}}} & \text{(6.4)} \\{A = ( {R^{2} + X^{2}} )^{1/2}} & \text{(6.5)} \\{\varphi = {{- \tan^{- 1}}\frac{X}{R}}} & \text{(6.6)}\end{matrix}$

[0207] [In Eqs (6.1) to (6.6), R and X indicate the real part and theimaginary part, respectively, and A and φ indicate the amplitude andphase delay, respectively.

[0208] R, X, A, and φ can be readily measured by a lock-in amplifier orthe like.

[0209] Since the relations (6.7) to (6.10) below hold, R, X, A, and φcan be written as in Eqs (6.11) to (6.14). $\begin{matrix}{\frac{\partial R_{i}}{\partial\omega} \equiv {- {\int_{0}^{\infty}{t_{i}{s(t)}{\exp \lbrack {- {\sum\limits_{i = 1}^{N}( {\mu_{ai}{ct}_{i}} )}} \rbrack}\sin \quad \omega \quad t{t}}}}} & \text{(6.7)} \\{\frac{\partial X_{i}}{\partial\omega} \equiv {- {\int_{0}^{\infty}{t_{i}{s(t)}{\exp \lbrack {- {\sum\limits_{i = 1}^{N}( {\mu_{ai}{ct}_{i}} )}} \rbrack}\cos \quad \omega \quad t{t}}}}} & \text{(6.8)} \\{\frac{{\partial\ln}\quad A_{i}}{\partial\omega} \equiv {{- \frac{1}{A^{2}}}( {{R\frac{\partial R_{i}}{\partial\omega}} + {X\frac{\partial X_{i}}{\partial\omega}}} )}} & \text{(6.9)} \\{\frac{\partial\varphi_{i}}{\partial\omega} \equiv {{- \frac{1}{A^{2}}}( {{R\frac{\partial X_{i}}{\partial\omega}} - {X\frac{\partial R_{i}}{\partial\omega}}} )}} & \text{(6.10)} \\{R = {{R( {\mu_{a} = 0} )} + {c{\sum\limits_{i = 1}^{N}{\int_{0}^{\mu_{ai}}{\frac{\partial X_{i}}{\partial\omega}{\mu_{a}}}}}}}} & \text{(6.11)}\end{matrix}$

[0210] The first term on the right side in each equation of Eqs (6.11)to (6.14) is a response on the assumption of no absorption, and thesecond term is a term dependent upon absorption.

[0211] As described above, there are four types of description methodsrepresented by Eqs (6.7) to (6.10) in the PMS and it is thus possible toimplement optical CT image reconstruction methods corresponding to therespective equations. On that occasion, their respective weightfunctions can be calculated using the values obtained by substitutingμ_(al)=μ_(ν) into Eqs (6.7) to (6.10) and using the relations of Eqs(6.11) to (6.14), in much the same manner as described previously.

[0212] The above detailed the preferred embodiments of the presentinvention, but it is noted that the present invention is by no meansintended to be limited to the above embodiments.

INDUSTRIAL APPLICABILITY

[0213] As described above, the present invention has permitted thedirect and quick calculation of the time-resolved path length li invoxel i by solving the forward problems for the imaginary medium assumedto have the same scattering characteristics and boundary conditions asthe measured medium and be nonabsorbing, so that the calculation timebecame shorter than in the conventional calculation of forward problemswith absorption. Further, the invention obviated the need for theconventional time-consuming forward problem calculation of varying theabsorption coefficient of each voxel and repeatedly solving the forwardproblems by the number of times equal to the total number of voxels.Accordingly, the present invention has permitted the provision of themethods capable of quickly calculating the photon path distributioninside the scattering medium of the heterogeneous system like livingtissue.

1. A method of calculating a photon path distribution inside ascattering medium, comprising: a first step of dividing a measuredmedium being the scattering medium, into N (1≦N) voxels each having apredetermined size and assumed to have a uniform absorption distributionin a measurement wavelength region; a second step of hypothesizing animaginary medium that can be assumed to have the same shape, boundaryconditions, and scattering distribution as the measured medium, have auniform refractive index distribution, and be nonabsorbing, for themeasured medium divided into the N voxels; a third step of setting alight injection position for injection of impulse light, and a lightdetection position for detection of an impulse response s(t) at a time tof the impulse light having been injected at the light injectionposition and having propagated in the imaginary medium, on a surface ofthe imaginary medium; a fourth step of calculating the impulse responses(t) detected at the light detection position after injection of theimpulse light at the light injection position into the imaginary medium;a fifth step of calculating a photon existence probability densityUi(t′,t) that an ensemble of photons constituting the impulse responses(t) detected at the time t at the light detection position afterinjection of the impulse light at the light injection position into theimaginary medium has existed in an arbitrary voxel i (i=1, 2, . . . , N;1≦N) in the imaginary medium at a time t′ (0≦t′≦t) before the detectionat the light detection position; a sixth step of calculating a photonexistence cumulative probability Ui(t) that the ensemble of photonsconstituting the impulse response s(t) detected at the time t at thelight detection position has existed in said voxel i, using the photonexistence probability density Ui(t′,t); and an eighth step ofcalculating a time-resolved path length li in the voxel i, using thephoton existence cumulative probability Ui(t) and the impulse responses(t).
 2. The method according to claim 1, further comprising, subsequentto the eighth step, a ninth step of calculating a weighted mean Li fortime of the time-resolved path length li calculated for each arbitraryvoxel i of the imaginary medium.
 3. The method according to claim 1,further comprising, subsequent to the eighth step, a tenth step ofcalculating a weight function Wi corresponding to a predeterminedmeasuring method for quantitatively measuring a concentrationdistribution of an absorptive constituent inside said measured medium,using said time-resolved path length li calculated for each arbitraryvoxel i of the imaginary medium.
 4. A method of calculating a photonpath distribution inside a scattering medium, comprising: a first stepof dividing a measured medium being the scattering medium, into N (1≦N)voxels each having a predetermined size and assumed to have a uniformabsorption distribution in a measurement wavelength region; a secondstep of hypothesizing an imaginary medium that can be assumed to havethe same shape, boundary conditions, and scattering distribution as themeasured medium, have a uniform refractive index distribution, and benonabsorbing, for the measured medium divided into the N voxels; a thirdstep of setting a light injection position for injection of impulselight, and a light detection position for detection of an impulseresponse s(t) at a time t of the impulse light having been injected atthe light injection position and having propagated in the imaginarymedium, on a surface of the imaginary medium; a fifth step ofcalculating a photon existence probability density Ui(t′,t) that anensemble of photons constituting the impulse response s(t) detected atthe time t at the light detection position after injection of theimpulse light at the light injection position into the imaginary mediumhas existed in an arbitrary voxel i (i=1, 2, . . . , N; 1≦N) in theimaginary medium at a time t′ (0≦t′≦t) before the detection at the lightdetection position; a sixth step of calculating a photon existencecumulative probability Ui(t) that the ensemble of photons constitutingthe impulse response s(t) detected at the time t at the light detectionposition has existed in said voxel i, using the photon existenceprobability density Ui(t′,t); a seventh step of calculating a sum U(t)of said photon existence cumulative probabilities Ui(t) for all thevoxels; and an eleventh step of calculating a time-resolved path lengthli in the voxel i, using the photon existence cumulative probabilityUi(t) and the sum U(t) of the photon existence cumulative probabilitiesUi(t) for all the voxels.
 5. The method according to claim 4, furthercomprising, subsequent to the eleventh step, a ninth step of calculatinga weighted mean Li for time of the time-resolved path length licalculated for each arbitrary voxel i of the imaginary medium.
 6. Themethod according to claim 4, further comprising, subsequent to theeleventh step, a tenth step of calculating a weight function Wicorresponding to a predetermined measuring method for quantitativelymeasuring a concentration distribution of an absorptive constituentinside said measured medium, using said time-resolved path length licalculated for each arbitrary voxel i of the imaginary medium.